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ABSTRACT 

We explore a systematic approach to the analysis of primordial non-Gaussianity using fluctuations in tem- 
perature and polarization of the Cosmic Microwave Background (CMB). Following Munshi & Heavens 
(2009), we define a set of power-spectra as compressed forms of the bispectrum and trispectrum derived 
from CMB temperature and polarization maps; these spectra compress the information content of the 
corresponding full multispectra and can be useful in constraining early Universe theories. We generalize 
the standard pseudo-C; estimators in such a way that they apply to these spectra involving both spin-0 
and spin-2 fields, developing explicit expressions which can be used in the practical implementation of 
these estimators. While these estimators are suboptimal, they are nevertheless unbiased and robust hence 
can provide useful diagnostic tests at a relatively small computational cost. We next consider approxi- 
mate inverse-covariance weighting of the data and construct a set of near-optimal estimators based on that 
approach. Instead of combining all available information from the entire set of mixed bi- or trispectra, 
i.e multispectra describing both temperature and polarization information, we provide analytical con- 
structions for individual estimators, associated with particular multispectra. The bias and scatter of these 
estimators can be computed using Monte-Carlo techniques. Finally, we provide estimators which are 
completely optimal for arbitrary scan strategies and involve inverse covariance weighting; we present the 
results of an error analysis performed using a Fisher-matrix formalism at both the one-point and two-point 
level. 
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INTRODUCTION 



The Cosmic Microwave Background (CMB) has the potential to provide cosmologists with the cleanest statistical characterization of primordial fluc- 
tuations. In most early universe studies the primordial fluctuations are assumed to be nearly Gaussian, but the q uest for an experime ntal detection of 
primordial non-Ga ussian is of considerable importance. Early observational work on the bispectrum from COBE l IKomatsu et alj|2002f) and MAXIMA 
( Santos et al. 2003) wa s followed by much more accurate analysis using data from the Wilkinson Microwav e Anisotropy Prob e (WMAPfl fkomatsu et all 
l2003l ; ICreminelli et alJl2007l ; ISpergel et alj2007l) .With the recent claim of a detection of non-Gaussianity jYadav & Wa ndelt 2008]) in the 5-year WMAP 
data, interest in non-Gaussianity has received a tremendous boost. However, these, and most other, studies of primordial non-Gaussianity focus primarily 
on data relating to temperature anisotropics whereas inclusion of E-polarization data could, in principle, increase the sensitivity of such tests. Ongoing 
surveys, such as that derived from the Planck Surveyofl will throw more light on this question both by improving sensitivity to both temperature and 
polarization fluctuations. It is the primary purpose of this paper to extend a number of previously obtained analytical results in order to design a new set 
of diagnostic tools for the analysis of primordial non-Gaussianity using using both temperature and polarization data. 

The CMB polarization field can be decomposed into a gradient part with even parity (the "electric" or E-mode) and a curl part with odd parity (the 
"magnetic" or B-mode), with important ramifications for its interpretation in a cosmological setting. Scalar (density) perturbations, at least in the linear 
regime, are unable to generate B-mode polarization directly. Secondary effects such as gravitational lensing can generate magnetic polarization from an 
initial purely electric type, but only on relatively small angular scales. Tensor (gravitational-wave) perturbations, however, can generate both iJ-mode and 
B-mode on large scales. Detection of £>-mode polarization on large angular scales is therefore, at least in principle, a tell-tale signature of the existence 
of gravitational waves predicted to be generated during inflation. The amplitude of these tensor perturbations also contains information relating to the 
energy scale of inflation and thus has considerable power to differentiate among alternative models of the early Universe. 



1 http://map.gsfc.nasa.gov/ 

2 http://www.rssd.esa. int/index.php?project=Planck 
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A number of experiments therefore are either being planned or currently underway to characterize the polarized CMB sky, including Planck. The 
primary statistical tools used currently to characterize the stochastic polarized component are the power spectra of the E and B-mode c ontributions; 
these have been studied in the literature in great detail for various survey strategies, non-uniform noise distribution, and partial sky coverage dHivon et al.l 
2002). Most studies in this vein involve a (computationally extensive) maximum likelihood analysis or a quadratic estimator, which have limited ability 
to handle maps with large numbers of pixels, or the so-called pseudo-C;s estimator (PCL) which uses a heuristic weighting of various pixels, instead of 

an optimal or near-optimal one. 

Information about non-Gaussianity, however, is not contained in the power-spectrum Cz, however it is estimated. A recent study by Mun shi et al.l 



(2009) extended the PCL approach to the study of non-Gaussianity by taking into account higher-order spectra, as well as clarifying the relationship of 
the estimators obtained to the optimal ones. This was done for temperature only, so one of the aims of this paper is to generalise this early work to take 
into account both temperature and polarization. 

Approaches based on PCL do not require detailed theoretical modelling of the target bispectrum. Optimal estimators on the other hand, work using 
a matched-filter approach which does require analytical m odelling. The simplest inflationary models - based on a single slowly-rolling scalar field - 
are associated w ith a very small level of non-Gaussianit y dSalopek & Bondll990b 1 199 lb iFalk et all ll 993b iGangui et all 1 1994 lAcouaviva et al.l 120031 ; 
Maldacena 2003); see lBartolo. M atarrese & Riottol d2006l) and references there in for more details. More elaborate variations on the inflationary theme, 



such as those involving multiple scala r fields i Lvth. Ungarelli & Wands 2003[). features in the inflationary potential, non-adiabatic fluctuations, non- 
standard kinetic terms, warm inflation dGupta. Berera & Heavensll2002b iMoss & Xiondl2007l) . or deviations from Bunch-Davies vacuum can all lead to 
much higher level of non-Gaussianity. 

Generally speaking, the apparatus used to model primordial non-Gaussianity ha s focused on a phenomen ological 'local /jvl' parametrization in 
terms of the perturbative non-linear coupling in the primordial curvature perturbation dKomatsu & S pergel 2 00lh : 

= <$> L (x) + f NL [$l(x) - (<f>l(x))}+g NL $l(x) + ..., (1) 

where <&l (x) denotes the linear Gaussian part of the Bardeen curvature ( Bartolo, Matarrese & Riottcl2006l) and /jvl is the non-linear coupling parameter. 
A number of models have non-Gaussianity which can be approximated by this form. The leading order non-Gaussianity therefore is at the level of the 
bispectrum or, equivalent, at the three-point level in configuration space. Many s tudies involving primordial non-G aussianity have used the bispectrum , 
motivated by the fact that it contains all the information about frvr, jB abich 2005). This has been extensively studied llKomatsu. Spergel & Wandeli2005b 
ICreminellj200l ; ICreminelli et al]2006l ; lMedeiros & Contaldoll2006l ; ICabella et al. l2006l ; lLiguori et al]2007blSmith. Senatore & Zaldarriagj2009l) . with 
most of these measurements providing convolved estimates of the bispectrum. Optimized 3-point estimators were introduced bylHeayensldl998l). and have 
been successively developed dKomatsu. Spergel & Wan delt 2005; ICreminelli et al.l2006l ; ICreminelli. Senatore. & Zaldarriagall2007blSmifh. Zahn & Dord 
2000; Smith & Zaldarriaga 2006) to the point where an est imator for Jnl which saturates the Cramer- Rao bound exists for partial sky coverage and 
inhomogeneous noise ( Sm ith. Senatore & Zaldar riaga 2009) . Approximate forms also exist for equilateral non-Gaus sianity, which may arise in models 
with non-minimal Lagrangian with higher-derivative terms dChen. Huang & Kachru 2006 ; Chen, Easth er & Liml2007h . In these models, the largest signal 
comes from spherical harmonic modes with l\ ~ £2 — £3, whereas for the local model, the signal is highest when one £ is much smaller than the other 
two. Moreover, the covariances associated with these power-spectra can be computed analytically for various models, thereby furnishing methods to test 
simulation pipeline in a relatively cost-effective way. 

As we mentioned above, most analysis of the CMB bispectrum take into account only temperature information because of lower sensitivity associated 
with polarisation measurements. However, the overall signal-to-noise ratio of a detection of non-Gaussianity can in principle be improved by incorporating 
.E-type polarisation information in a joint analysis. Several ground-based or future experiments have either started or will be measuring _E-polarisation, 
and ongoing all-sky experiments such as Planck w ill certainly improve the signal to noise. Moreover, there are many planned experiments which will 
survey the sky with even better sensitivity; see e.g. dBaumann et al.|[2009l) for the CMBPol mission concept study. It is therefore clearly timely to update 
and upgrade the available estimators which can analyze non-Gaussianity including both temperature and polarization. The main question is to how to do 
it optimally. 

The fast estimator introduced initially by iKomatsu, Spergel & WandelJ (2005) can handle temperature and polarization data. Extension of this 
estima tor to take into account a linear term, which can reduce the variance in the presence of partial sky coverage, was introduced later in lCreminelli et al] 
(2006). In the absence of a linear term the scatte r associated with the estima tor increases at high -£. A more general approach, that includes inverse 
variance weighting of the data, was introduced by ISmith"& Zaldarriag3 d2006l) and can make the estimator optimal. However most of these estimators 
tend to compress all available information into a single point. This has the advantage of reducing the scatter in the estimator but it throws away a great 
deal of potent ially relevant information. 

Recently iMunshi & Heavens! d20ld) have extended this approach by devising a method which can map the bispectrum, or even higher-order spec- 
tra such as the trispectrum, to associated power-spectra; these power s pectra are related to the cumulant correlators used in the context of projected 
galaxy surveys (Munshi, Melott & Coles 2000; Szapudi & Szalav 1999). These compressed spectra do not contain all the information contained in the 
multispectra but they certainly do carry more information than a single number, and for certain purposes, such as /ml estimation, can be completely 
optimal. 

In this paper we extend the analysis o f IMunshi & Heaven j ( 120101) to the case of a joint analysis of temperature and iJ-mode polarization data. We 
provide realistic estimators for power spectra related to mixed bispectra and trispectra (i.e. multispectra incorporating and describing both temperature and 
polarization properties of the radiation field) in the presence of partial sky coverage. We then generalize these results to the case of optimized estimators 
which can handle all realistic complications by an appropriate optimal weighting of the data. We also provide a detailed analysis of pseudo-C« based 
estimators. Though suboptimal, these estimators are unbiased and turn out to be valuable diagnostics for analyzing large number of simulations very 
quickly. 

The plan of the paper is as follows. In §2 we describe the basic results needed for the PCL analysis. We then introduce the power spectrum associated 



3 See e.g. http://cmbpol.uchicago.edu/workshops/path2009/abstracts.html for various ongoing and planned surveys. 
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with the mixed bispectrum. The formalism is sufficiently general to handle both E and £?-mode, but for simplicity and practical relevance we give result s 
for the _E-mode only. These results are next generalized in §4 to near-optimal estimators which were introduced bv lKomatsu. Spergel & WandelJ d2005t) . 
The optimization depends on matched-filtering techniques based on theoretical modeling of primordial non-Gaussianity presented in §3 and relies on 
Monte-Carlo standardization of bias and scatter. Next, we consider the approach where various fields are weighted by inverse covariance matrices to 
make them optimal in §5. Though this approach does not rely on Monte-Carlo simulations, accurate modelling of inverse covariance matrix can be 
computationally expensive and can only be achieved for maps with relatively low resolution. 



2 METHOD OF DIRECT INVERSION OR PSEUDO C' e APPROACH 

Maximum likelihood analysis techniques, or related quadratic estimators, are generally used to estimate the power spectrum from a given cosmological 
data set such as a map of the CMB. This approach relies on (optimal) inverse-covariance weighting of the data which, unfortunately, is virtually impossible 
to implement in practice for large data sets. However, faster direct i nversion techniqu es with heuristic weighting schemes are faster and are typically 
employed for estimation of power spectrum from cosmological data dHivon et al. 2002). The alternative, PCL, technique discussed above can, however, 
be used to study the power spectrum related to the bispectrum. This method is not optimal, but remains unbiased and can act as a precursor to more 
computationally expensive studies using various optimal estimators. 

This section is devoted to the generalization of the pseudo-C; (PCL) results to the skew-spectrum associated with spinorial fields such as polarization 
radiation. The results we shall derive are valid for all-sky coverage as well as for partial sky coverage; we relate the two using a coupling matrix which 
depends on the power spectrum associated the mask describing the coverage. We provide results for the coupling matrix for various generic combinations 
of spin-0 and spin-2 fields. These results extend the results obtained previously for the temperature-only case (spin-0). 



2.1 The Power Spectrum associated with the Bispectrum 

We start by defining the complex field P± (Cl) constructed from the Stoke's parameters Q(Cl) and U (Cl). 

P±(Q,) = Q(Cl)+iU(Cl) (2) 

The fields denoted P± (Cl) are spin-2 (tensor) fields, whereas the temperature field can be represented as a spin-0 (scalar) field on the surface of the sky. 
The appropriate multipole expansion of these objects is performed using the spin-2 spherical harmonics, ±2Yi m , as basis functions. We will denote the 
spin harmonics of spin s with s Yi m (Ci) on the surface of a unit sphere. 

^2Yi m (Q.)(Ei m ± Bi m ); [P±]im = {Eim ± Bi m )- (3) 

lm 

The terms Ei m and Bi m are the harmonic components of Electric and Magnetic components respectively. It should be clear that the fields P±, constructed 
from Q and U, correspond to spin =p2 respectively. 

Estimation of power spectrum for the E and B fields from experiments with partial sky coverage is of great importance for cosmological experiments 
and has attracted a great deal of interest. 

We start by considering two fields on the surface of a sphere X(Cl) and y(Cl), which are respectively of spin x and y. These objects can be P± 
or ST. Resulting product fields such as P+P- can be of spin zero, while those like P+ and P* are of spin —4 and +4 respectively. The results can 
also involve fields such as 8tP+ which is a spin —2 field. All these spin-s fields can be decomposed using spin-harmonics a Yi m (Cl) given above. These 
spin harmonics are generalization of ordinary spherical harmonics Yj m (Cl) which are used to decompose the spin-0 (or scalar) functions e.g. 5T defined 
over a surface of the celestial sphere. Throughout we will be using lower case symbols x, y to denote the spins associated with the corresponding fields 
denoted in italics. The fields that are constructed from various powers of P+ and P_ can be expanded in terms of the corresponding spin-harmonics basis 
s Yi m (Cl). The product field such as [^(0)^(0)] which is of spin x + y can therefore be expanded in terms of the harmonics x + y Yi m (Cl). The following 
relationship therefore expresses the harmonics of the product field in terms of the harmonics of individual fields. 



[x(fi)y(fi)] lm = J dsix(si)y(ti)[ x+y Y lm ({i)] (4) 
^^ imi ^ 2m2 y , ty iimi (n)][ a y i2m2 (n)][, +M y^(n)]dn; x,y&p ± (5) 



I h k l Ii h l V . , / (2l 1 + l)(2Z 2 + l)(2Z + l) 

2^x hmi y hm2 i lxhl \^ x y _ {x + y) )^ mi ma _ m j, /j lW = y ^ • (6) 



The expression derived above is valid for all-sky coverage; it will be generalized later to take into account arbitrary partial sky coverage. We have assumed 
that any noise contamination is Gaussian, so that it will not contribute to this non-Gaussianity statistic. To define the associated power spectra we can 
write: 



C i -2mL[ ly W»--L B 'i'W (2T+T) (x y -(x + y) ■ (7) 

Here we have introduced the bispectrum B^ff which can be related to the harmonics of the relevant fields by the following equation: 
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B? y f — N (Xi irn ,yi 2m2 Zi 3m „) c ( 3 ) . (8) 



The matrices denote the Wigner-3j symbols dEdmondslll9"6il which are on ly non-zero when the quantum numbers U and rm satisfy c e rtain c on 



Cooravl < f200lh valid for the case of temperature; see also lChen & Szapudi f2007h for 



ditions. For x = y = these results generalize those obtained in 
related discussions. The result derived above is valid for a general E and B type polarization field. The contribution from B-type magnetic polarization 
is believed to be considerably smaller than the i5-type electric polarization. The results derived above will simplify considerably if we ignore the B type 
polarisation field in our analysis. In general the power-spectra described above will be complex functions, though the real and imaginary parts can be 
separated by considering different components of the bispectrum. However, if we assume that the magnetic part of the polarization is zero the following 
equalities will hold: 

tjSTSTST _ r>TTT V) P ± STST _ ^ETT I3 P±P±ST _ jjEET r> p ± p ± p ± _ r>EEE tQ . 

The result detailed above is valid for all-sky coverage. Clearly, for our results to be relevant in practical applications, we need to add a Galactic mask. 
If we consider the masked harmonics associated where the (arbitrary) mask w(Cl) then, expanding the product field in the presence of the mask, we can 
write: 



[x(Ci)y(Ci)w(Ci)]i m = ^2 x h mi yi 2m2 wi ama /[»y !l m 1 (ft)]| 1( yj a m a (A)][i r i B T^][ 1 .-^yj^]<iA 

= J2J2^ V^»U".Ui ( t y _(J+y))(m 1 m 2 - m > ) J l-i*+v) Y i'm>]Yi ama [*+vYi* m ] 

= — Xhm i yi 2 -m 2 'Wlam a Il 1 l 2 l'Il<la,l 



h h V \ h h V 
x y ~(x + y) I [ mi ffi2 — m' 

h , .., n J, ... 1 ( L h L V do) 



{x + y) — (x + y) I \ m' m a — m 

In simplifying the relations derived in this section we have used the relationship Eq.d34t and Eq.J35b. We have also used the fact that S Y^ = 
( — l) m+s Yi,~ m , where * denotes the complex conjugate. Similarly, we can express the pseudo-harmonics of the field Z(Q) observed with the same 
mask in terms of its all-sky harmonics. 



[Z(Cl)w(Cl)] lm = / dCl[Z(Cl)w(Cl)][ z Y; m ] = Z hm3 w hmb J [ z Y hm3 (Cl)}[Y lbmb (n)}[ z Yr m }dCl 

liTTli 

= 12 53^3-3 W '»m b Il 3 l b l ( ^ ^ l _ z J f ^ ^ J m J • (ID 

The harmonics of the composite field [A"(f2)3^(n)] when constructed on a partial sky are also functions of the harmonics of the mask used, wi m - 
The simplest example of a mask would be w — 1 within the observed part of the sky and w = outside. For a more complicated mask, the harmonics 
wim are constructed out of spherical harmonics transforms. We also need to apply the mask to the third field which we will be using in our construction 
of the power spectrum related to the bispectrum (which we sometimes refer to as the skew-spectrum) associated with these three fields. We will take the 
masks in each case to be the same, but the results could be very easily generalized for two different masks. The pseudo power spectrum is constructed 
from the masked harmonics of the relevant fields in the form of a cross-correlation power spectrum: 



c? y > z = [x(Ci)y(Ci)w(Ci)] lm [z(Ci)w(Ci)]t m . (12) 

m— — I 

The pseudo power spectrum Ci is thus a linear combination of the modes of its all-sky counterpart Ci. In this sense, masking introduces a coupling of 
modes of various order which is absent in the case of all-sky coverage. The matrix M H i encodes the information regarding the mode-mode coupling and 
depends on the power spectrum of the mask wi. For surveys with partial sky coverage the matrix is not invertible which signifies the loss of information 
due to masking. Therefore, binning of the pseudo skew-spectrum or any higher-order version based on, e.g., kurtosis, may be necessary before it can be 
inverted, which leads to the recovery of (unbiased) binned spectrum. 

&i y ' Z = M w' ZC v y ' Z - (°) 
V 

The all-sky power spectrum C* y ' z can be recovered by inverting the equation related to the accompanying bispectrum by Eq.(j7j discussed above. The 
mode-mode coupling matrix depends not only on the power spectrum |ioj| of the mask but also on the spin associated with the various fields which are 
being probed in construction of skew-spectrum: 
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^4E(tf + D(a. + i)(,i B l o _ ( /+v))G o -O^l 2 (14) 

This expression reduces to that of the temperature bispectrum if we set all the spins to be zero x = y — z = 0, in which case the coupling of various 
modes due to partial sky coverage only depends on the power spectrum of the mask. In case of temperature-only (spin-0) analysis we have seen that the 
mode-mode coupling matrix do not depend on the order of the statistics. The skew-spectrum or the power spectrum related to bispectrum, as well as its 
higher order counterparts, such as the power spectrum related to trispectrum, can all have the same mode-mode coupling matrix in the presence of partial 
sky coverage. However, this is not the case for power-spectra related to the polarization multispectra; that depends on the spin of various fields used to 
construct the bispectrum. It is customary to define a single number associated with each of these bispectra. The skewness is a weighted sum of the power 
spectrum related to the bispectrum S* y ' z = y^,{2l + l)C^ y ' z . 

We list the specific cases of interest below. The rela tions between the bispectra and the associated power spectra generalize cases previously obtained 
where only the temperature bispectrum was considered (Coorav 2001). 



{2k + 1){21 2 + 1) k h I \ n „ 
C i =L B Wt/ ^{21 + 1) { ) (15) 

M^ = iJ2^' + mu + i) ( l <o 5 ) ( \ \ U ) kj 2 (i6) 

s^yT* E E T E E 

Next we can express the power spectrum C\ ' probing the mixed bispectrum Bi t t by the following relation: 



hh 



r TE,E _V- WTEE / (2i 1 +l)(2i 2 + l) f k h I \ nr . 

^ -L^'y 4n{2l + l) I 2 -2 I (17) 



M^ = lJ2{2i' + i){2i a + i)( I l ; M(J l ; <;W. as) 



E E E E E E 

Similarly, for the case C l ' , which probes the bispectrum B lll2l , we have the following expressions. 



n EE,E _ 'ST^ nEEE /(2Zi + l)(2Z 2 + l) k h I \ nm 

C ' ' 4n{2l + l) ^4 -4 ) (19) 

«T = iD«' + 1 K !a - + 1 )(l o - 4 ) (2 -2)^ I 2 - ^ 

Other power spectra, such as Cf E ' T and those involving B-mode polarization, can be expressed in terms of the underlying bispectrum in a very 
similar manner. It is likely that future high sensitivity experiments can probe these power spectra individually. As these are unbiased estimators they 
can also be useful in testing simulation pipelines involving a large number of Monte-Carlo realizations, which are routinely used for standardization and 
characterization of data-analysis pipelines. The inversion of the coupling matrix M H i can require binning, in which case the relevant binned coupling 
matrix Mw in the case of small sky coverage (such as ground-based and balloon-borne) surveys will lead to recovery of binned power spectra Ci b . The 
flat sky analogs of these results can be useful in other cosmological studies, including weak lensing observations, involving spin-2 fields (mapping out 
the shear or 7 on the surface of the sky) that covers a small fraction of the sky. We plan to present results of such an analysis elsewhere (Munshi et al. 
2010, in preparation). 



2.2 The Power spectrum associated with the Trispecrum 

The bispectrum represents the lowest-order deviation from Gaussianity; the next highest order is the trispectrum. There are many reasons for wanting to 
go beyond the lowest-order description. Many studies relating to secondary anisotropies have shown that ongoing surveys such as Planck can provide 
information about non-Gaussianity beyond lowest order. These include the mode coupling of CMB du e to weak lensing of CMB as well as the other 
secondary effects such as thermal Sunyaev-Zeldovich (tSZ) and kinetic Sunyaev-Zeldovich effect (kSZ) JCoorav & Hu 2000). With the recent (claimed) 
detection of primordial non-Gaussianity in the CMB, there has also been a renewed interest in detection of primordial non-Gaussianity beyond lowest order 
with reasonable signal-to-noise. Such studies will constrain the parameters gNL and /ml independently. The signal-to-noise for such constraints is weaker 
compared to the ones achieved using bispectrum, but with anticipated future increases in sensitivity of ongoing as well as fut ure planned surveys such 
studie s will play an import ant role in fu t ure. F uture 21cm surveys can also provide valuable information about the trispectrum dCoorav. Li & Melchiorril 
120081) . A recent study bv lMunshi et~al] J2009h constructed an estimator for trispectrum which c an work with te mperature data and which is optimized 
to detect primordial non-Gaussianity. This statistic was also applied to WMAP 5-year data in JSmidt et"al]|2010]) to provide the first constraints on tnl 
and <7nl, the third order corrections to primordial perturbations in non-Gaussian models. We generalise and extend these works here to include both 
temperature and polarizations data for independent or joint estimate of such quantities. 
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2.2.1 The Two-to-Two Power spectrum 

In constructing the power-spectrum related to the trispectrum, we start with two fields U and V respectively on the surface of the sphere. As before, we 
can take specific examples where these fields are either P± or ST. We will keep the analysis generic here, but will consider more specific examples later 
on. The spins associated with various fields are denoted by lower case symbols, i.e. u and v. The product field now can be expanded in terms of the spin 
harmonics of spin u + v, as was done in Eq.©. Similarly, decomposing the other set of product field we obtain [W(n)<%"(n)]i m . We next construct the 
the power spectrum associated with the trispectrum from these harmonics: 



r <uv,wx _ 
°' -21 



— [K(Ci)V(Ci)]i m [W(Ci)x(Ci)]t m . (2i) 



This particular type of p ower spectrum associated with trispectra has been studied extensively in the literature in the context of CMB studies (see e.g. 
ICoorav & Kesdenl d2003h ). It is one of the two degenerate power spectra associated with trispectrum. After going through very similar algebra outlined in 
the previous section we can express the C l v ' in terms of the relevant trispectra which it is probing. The resulting expression, in the absence of any 
mask, takes the following form: 



^1 ^2 >H 



r uv.wx _ ST T U >i v h m ( h h 1 \ f la li 1 \ / (2ii + l) (2k + l) / (2b + l)(2i4 + l) nj . 

~ 2^ J w,,x u Wl u v ~( u + v ) )\ x y -(x + y) )\j (21 + 1) V (21 + 1) ' K ' 



U,V, W,X 6 P±,St; 



LM v 7 v 7 

In the presence of a completely general mask w(Q) the pseudo-C^s or PCLs will have to be modified. This involves computing the spherical harmonics 
of the masked field U(Cl)V(Cl)w(Cl) and cross-correlating it against the harmonics of W(Cl)X(Cl)w(Cl). 



c: 



i i 

■uv,wx _ 1 \ ^ 



21 + 1 ^ 

m— — I 



[UV w] lm [WX w]* lm ; wi = — - E wi m w* lm . (24) 



Again, the resulting PCLs are linear combinations of their all-sky counterparts. The mixing matrix which encodes information about the mode 
mixing will depend on the power spectrum of the mask as well as the spins of all four associated fields. The mixing matrix M U i expressed in terms of the 
Wigner's 3j symbols takes the following form: 



M uv,*y = 1 \^( 2 f' + l)(2f + l) f * l ° /',)([ n t 1 '. ^ ) u,v,x,y e ±2,0. (25) 

47r ' y u + v —(u + v) J \ x + y — (x + y) 1 



The pseudo-C^ s expressed as a linear combination of all-sky power spectra can now be expressed using the following relationship: 



d uv,wx = y M u r , xc uv,wx^ (26) 



For nearly complete sky surveys, and with proper binning, the mixing matrix Mu< can be made invertible. This provides a unique way to estimate 
all-sky C,, ' and extract the information it contains about the trispectra. For a given theoretical prediction, the all-sky power spectra C^, vw ' x can 
be analytically computed. Knowing the detai led model of an exp erimental mask then allows us to compute the observed C, ' accurately. The results 
presented here generalizes those obtained in dMunshi et aI1l2009l) for the case of polarisation studies. 

These results assumes generic field variables X(Q), y(Cl) which can have arbitrary spin associated to them. We next specify certain specific cases 
where we identify three of the fields X, y, Z = St and Z = P+ — E. The other combinations can also be obtained in a similar manner. 



r TT,TE _ ST T Tl ^(i\ h h l\ h h I \ / (2Ii + l)(2fa + l) / (2i 3 + l)(2i 4 + l) , 97 , 

° l Vu (,) l o o o 2 o -2 jy (2/TI) V (2T+1) ' 

^1 >^2^3>^4 



Mr 2 = ^E^' + 1 )^ + 1 )( I 2 -2 )'^' 2 = 



Other estimators for mixed trispectra involving different combinations of B-polarization and temperature anisotropy St can be derived in a similar manner 
and can provide independent information of corresponding trispectra. As before, we have ignored the presence of B-type polarization in our analysis. 
The presence of a non-zero B-mode can be dealt with very easily in our framework but the resulting expressions will be more complicated. 



2.2.2 The Three-to-One Power spectrum 

The number of power spectra that can be associated with a given multispectrum depends on number of different ways the order of the multispectrum can 
be decomposed into a pair of integers. The bispectrum being of order three can be decomposed uniquely 3 = 2+1 and has only one associated power 
spectrum. On the other hand the order of the trispectrum can be decomposed in a two different ways, i.e. 4 = 3 + 1 = 2 + 2. Hence the trispectrum of 
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Figure 1. a(r) (lower panel) and /3(r) (upper panel) for Temperature are plotted as a function of I. Plots are based on WMAP7 parameters lLarson et alj201Ch 




Figure 2. Same as previous plot but for E-type polarisation. 



a specific type generates a pair of two different power spectrum associated with it; see, e.g. iMunshi et all d2009h for more details. The results presented 
here are generalizations for the case of non-zero spins. 

The other power spectrum associated with the trispectrum is constructed by cross-correlating product of three different fields [WVW] with the 
remaining field [X]. The cross-correlation power spectrum in terms of the multipoles is given by the following expression: 



[UVW]i 



[u(n)v(n)x(Ci)][ u+v+w Y lm (ny]dn- 



[x(n)}[ z Y tm (Q,y}dn; c; 



■uvw,x 



— [UVW]i m {X]t m . (29) 



; 21 

m—~l 

By repeated use of the expressions Eq.<|341> or, equivalently, Eq.d35t to simplify the harmonics of the product field in terms of the individual harmonics, 
we can express the all-sky result in the following form: 



h h L 
u v — (u + v) 



L l s I' 
(it + v) w — (u + V + w) 



(2L + 1) 



(2h + 1)(2Z 2 + 1) /(2L + l)(2/ 3 + l) 



(21 + 1) 



(30) 
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In the case of of partial sky coverage with a generic mask w(Cl) the relevant expression for the Pseudo-Cf s will as usual involve the harmonic transform 
of the mask. 



[uvww] lm = J [w(n)v(«)W(n)«;(n)][» +t)+w y4(fi)]dn; [x w ] lm = J [x(si)w(Q)]\„Y{^(Q)]dSi; =M u ,c% vm -* (3i) 

The mixing matrix has the following expression in terms of various spins involved and the power spectra of the mask introduced above. Note that the 
mixing matrix is of different form compared to what we obtained for two-to-one power spectra. This is related to how various fields with different spins 
were combined to construct these two estimators. In case of temperature trispectrum (spin-0) the two mixing matrices take the same form. 

M u r ,* = 1 Y(2l' + l)(2l a + l)( . f / , A( 1 IK I 2 ; (32) 

11 4tt y (u + v + w) — (u + v + w) J \ x —x J 

The expressions derived here apply to a general mask and to correlated Gaussian instrument noise; any non-Gaussianity from the noise will have to 
be allowed for. These calculations lead to optimal estimators, but the weights required depend on the precise model of non-Gaussianity being assumed. 
We present the models for primordial non-Gaussianity in the next section and use them to construct optimal estimator in our later discussions. 

For a specific example we choose U = V = W = St and X — £ . In this case, the three-to-one estimator takes the following form: 



c ttt,e = <y T T ' lT "(L)f h h L \( L h 1 \ / (2ii + l)(2k + l) / (2L + l)(2i 3 + l) - 



T l3 E t \ / i o I \ I Y (2L + 1) V (21 + 1) 



M^ = iY,w+mi«+v ( 5 u l ; U ) i^i 2 - ^ 

The bispectrum is defined through a triangular configuration in the multipole space, the trispectrum with a quadrilateral (which can be decomposed 
into two constituent triangles). The Wigner-3j symbols enforce these constraints at various levels for bispectrum and trispectrum. The two-to-one power 
spectrum or skew-spectrum introduced in previous sections essentially sums over all possible configurations of the triangle obtained by keeping one of 
its sides fixed. In an analogous manner the two probes of trispectrum introduced here are linked with two different configuration in the harmonic domain. 
The two-to-two power spectrum keeps the diagonal of the quadrilateral fixed whereas the three-to-one power spectra keeps one of the sides fixed while 
summing over all other possible configurations. 

It is possible to introduce a window optimized to search selectively for information for a specific configuration for either the bispectrum or the 
trispectrum. However such analysis which introduces mode-mode coupling can not be generalized to the arbitrary partial sky coverage as there is already 
a coupling of modes because of non-uniform coverage of the sky. 

Unlike the bispectrum, the trispectrum has a non- vanishing contribution from the Gaussian component of a CMB map (in the same way that, while 
the third moment of a univariate Gaussian about the mean is zero, the corresponding fourth-order moment is not). The Gaussian contribution represents 
the unconnected component of the total trispectrum which needs to be subtracted out. This can be done by generating Monte-Carlo maps using an identical 
mask. The resulting Gaussian part of the sp ectra can then be su btracted from the estimates from the real data. The procedure is similar to the analysis of 
temperature-only data; for more details see iMunshi et al1d2009h . Noise (assumed Gaussian) can also be subtracted out following a similar procedure. The 
treatment requires a hit-count map from a realistic scanning strategy, for non-uniform distribution of noise, from where the variance of noise distribution 
in each pixel can be constructed. 

To derive and simplify the above expressions we have used the following results related to the overlap integral involving three spin harmonics which 
generalizes the well-known Gaunt Integrals involving spin-0 spherical harmonics. 



v «»■ v o \ <> „> ; (2Z + l)(2l' + l)(2l" + l) f I I' I" \ I I 1 I" \ 

,Y lm {n) s ,Y l , m ,{n) s „Y l „ m ,,(n)<m = x j — I m , m „ I ( s _ s > _ s » )■ (34) 

The above result can be cast in following form which is useful for expressing integral of more than three spin spherical harmonics in terms of Wigner-3 j 
symbols: 



v (6\ v (6\ ST v* ,/ (2£ + l)(2i' + l)(2i" + l) f L V I" \ L V I" , 

s Y lm (Cl) s ,Y lm (Cl) - ^ sY LM J I M m , m „ II _ s ). (35) 

LSM \ / \ 



-s —s 



The Wigner-3j symbols satisfy an orthogonality relationship which can be written as: 

I h h la \ ( h h h \ Si a i b S mamb 



£-~t \ mi m2 m a J \ mi 7712 m b ) 2l a + 1 

m\m2 



(36) 



In the remainder of the paper we will generalize the estimators by including inverse variance weighting. We will first show how to include the effect 
of noise and mask by using Monte-Carlo simulations. Next we will consider the exact inverse variance weighting where the inverse covariance matrix 
incorporates the effect of noise and mask. 

Note that incorporating a non-zero B-type (magnetic) polarisation is possible in our formalism presented above by conceptually trivial extension. 
We have only taken into account E-type (electric) polarisation to keep the derivations simple as E-modes are assumed to dominate presence of B-mode 
polarisation. 
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3 THE BISPECTRUM AND TRISPECTRUM WITH PRIMORDIAL NON-GAUSSIANITY 

The optimization techniques that we will introduce in this section follow the discussion in iMunshi & Heavens! J2010h and lMunshi etafl d2009h . The 
optimization procedure depends on construction of three-dimensional fields from the harm onic components of the temperatu re fields a; m with suitable 
weighting with respective functions (a, J3, /i) which describes primordial non-Gaussianitv dYadav. Komatsu. & Wandelt 200^). These weights make the 
estimators act in an optimal way and the matched filtering technique adopted ensures that the response to the observed non-Gaussianity is maximum when 
it matches with the target primordial non-Gaussianity corresponding to the weights. 

In the linear regime, the curvature perturbations which generate the fluctuations in the CMB sky are written as: 

af m =4n(-i) l J ^L$(k)A?(k)Y lm (ky, XeT,E. (37) 
We will need the following functions to construct the harmonic space trispectrum as well as to generate weights for the construction of optimal estimators: 

<**{?) = - / k 2 dkA? (fc)ji(fcr); A*(r) = - / k 2 dkP^{k)A? ! (*)ji(*r)i (r) = - / k 2 dk&? \k)j,(kr) (38) 
71 Jo * Jo 71 Jo 

For a more complete description of the method for extracting the trispectrum from a given inflationary model, see iHu] fcOOChlHu & Okamotol d2002h . 
In the limit of low multipoles, where the perturbations are mainly dominated by Sachs-Wolfe Effect the transfer functions A; (k) take a rather simple 
form A; (fc) = ji (kr, )/3 where r* = (rjo — rjdec) denotes the time elapsed between cosmic recombination and the present epoch. In general the transfer 
function needs to be computed numerically. The local model for the primordial bispectrum and trispectrum can be constructed by going beyond linear 
theory in the expansion of the Additional parameters f^L and ijjvl are introduced which need to be estimated from observation. As discussed in 

the introduction, qnl can be linked to r the scalar to tensor ratio in a specific inflationary model and hence expected to be small. 

$(x) = $ L (x) + f NL ($l(x) - <$!(x)» + <?iVL<£-!(x) + h NL (#t(x) - 3{$!(x)» + . . . (39) 

We will only consider the local model of primordial non-Gaussianity in this paper and only adiabatic initial perturbations. More complicated cases of 
primordial non-Ga ussianity will be dealt with elsewhere. In terms of the inflationary potential V(4>) associated with a scalar potential (f> one can express 
these constants as jHull200Cl) : 



5 1 d 2 \nV{<t>) 25 

/NL = 7777 7775 ' 3NL = 



6 8-kG d<j> 2 ' 54(8ttG) 2 



d 2 In V{(j>) ^ _ ^ d 3 \nV(4>) ^ ( dlnV(<t>) \ 



(40) 



There are two contributions to the primordial non-Gaussianity. The first part is parametrized by /nl and the second contribution is proportional to a 
new paramet er which appears at fourth or der which we de note b y <7nl- From theoretical considerations in generic models of inflation one would expect 
<?nl < r/50 JSeerv. L idsev & Sloth 2008). Following lHvJ |2000) we can expand the above expression in Fourier space to write: 

$2(fe) = / ^ <E, 4k + k 1 )$i(k 1 )-(2^) 3 fo(k) J ^Lp M ( kl ); $ 3 (fc) = J 7 ^3$4ki)<i' ; (k 2 )<E.r(k 1 ). (41) 

The resulting trispectrum T$ for the potential $ associated with these perturbations can be expressed as: 

/rf 3 K 
— - 5d (ki + k 2 - K)S D (ks + k 4 + K) T* (ki , k 2 , k 3 , k 4 , K) ; (42) 

where the reduced trispectrum T(ki, k 2 , k 3 , k 4 , K) can be decomposed into two distinct contributions: 

ri 2) (ki,k 2 ,k 3 ,k 4 ,K) = 4/£ L P (K)P*(ki)P*(k 3 ); Ti 3, (ki,k 2 ,k 3 ,k 4 ,K) = g NL {P (K)P^(ki)P4,(k 3 ) + P (K)P (ki)P*(k 3 ).} (43) 
The resulting mixed trispectrum, involving temperature and P-type polarization now can be written as: 



T yVi 3 X^( L ) — AfxthlhlaLlhliL J r 2 dnr 2 dr2F L (ri,r2)a{ l 1 (ri)l3i 2 (ri)a^(r2)l3^ 4 (r2) 

+gNX J I hhL Ii 3 i i L J r 2 dr^(r)<(r)[^(r)/?^(r) + ^(r)^(r)]. (44) 

For detailed descriptions of objects such as a, /3, fx, F L see fau & Okamotdl2002l : lH^l200d : lKorriatsu & Spergell200ll:lKogo et al.l2006l). The CMB 
bispectrum which desc ribes departures fr om Gaussianity at the lowest level can be analysed in a similar fashion; see IMunshi & Heavens! j2010r) for a 
detailed discussion and lSmidt et alj J2009h for a measurement in data. The corresponding expression is: 

B hlX = 2/nl/ W3 / r 2 dr [aft {r)tf 2 (r)P^(r) + a% (r)A^(r)^ (r) + a^(r)^ (r)# (r)] . (45) 

The extension to orders beyond those presented here involves higher-order Taylor coefficients and may not be practically useful as detector noise and 
the cosmic variance may prohibit any reasonable signal-to-noise. Numerical evaluations of the functions a, ft, fi can be performed by using the publicly 
available Boltzmann solvers such as C AME0 or CMBFAST0. 



4 http://camb.info 

5 http://www.cmbfast.org 
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Figure 3. Various bispectrum -related power spectrum (21 + 1)/Q 2 ' ^ plotted as function of angular scale I. We use /jvl = 1 for each of these plots. Plots are based on 
WMAP7 jLarson et alj20ld) parameters equated out to lmax = 500. 



A Gaussian fluctuation field has zero bispectrum. However, even a Gaussian map will have non-zero trispectrum; this corresponds to the unconnected 
part of the trispectrum. Its contribution can be expressed in terms of the cross power spectra associated with contributing fields: 

G% V * (L) = (-l) (,1+ ' 3 V( 2 *i + 1)(2Z 3 + l)Cli v C% x 5 hl2 S hu 6 L0 + (2L + 1) \{-l) h+h+L C?™ C?* 5 hh 5 l2li + C?* C? 2 W 5 h iji 2 i 3 ] (46) 

In the case of Gaussian maps, the trispectrum can be expressed completely in terms of the relevant cross power -spectra Cj of corresponding fields 
U and V. The estimators designed in later sections estimate the combined skew- or kurt-spectra and the Gaussian contributions are subtracted accordingly. 
The Gaussian maps that are used for Monte-Carlo estimates of bias and scatter are constructed to have same power spectrum as the non-Gaussian maps 
being analysed. 



4 PARTIAL SKY COVERAGE AND INHOMOGENEOUS NOISE 

In this section, we consider the inverse-variance weighting of the data. We will consider the all-sky case first and then introduce the analytical results 
that can handle data in the presence of partial sky coverage as well as correlated Gaussian noise. The method developed here relies on Monte-Carlo 
simulations to model observational artefacts. These estimators uses a weighted version of square temperature - temperature (or two-to-one) angular power 
spectrum. In a manner similar to its simpler version introduced in the previous section, this power spectrum extracts information from the bispectrum as 
a function of length of one side of the triangle in harmonic space, while summing over all possible configuration given by the change of other two sides 
of the triangle. 



4.1 Bispectrum-related Power Spectrum or Skew-spectrum 

Following iKomatsu. Spergel & WandelJ J2005t) . we first construct the 3D fields A(r, Cl) and B(r, Cl) from the expansion coefficients of the observed 
CMB map, aim- The harmonics here Ai m (r) and Bi m (r) are simply weighted spherical harmonics of the temperature field ai m with weights constructed 
from the CMB power spectrum C\ and the functions aj(r) and f3i(r) respectively: 

A u (r,Cl) =^Y lrn (£l)AY m (r); Af m (r) = aftr) ^^[C-^Mjw = c^VL; (47) 

lm hi' I'm-' 

B u {r,Cl) = J2Yi m (Cl)BY m (r); Bf m {r) = ff(r) ^ ^[CT^MjC = PYbi&Y m . (48) 

lm hi' I'm' 

The function bi represents beam smoothing, and from here onward we will absorb it into the harmonic transforms. The matrix C depends on both 
temperature (T) and E-type polarization power spectra Cj T and Cf , The cross-correlation power spectrum is denoted by C; X . 

cr k E ) ; [c]rl = ~k [ cr -k T ) ; ~ aU = [clwi,a " • (49) 

The determinant T>i introdu ced above is a function of the relevan t three power spectra introduced above T>i — Cj T Cf E — (C; ) 2 for joint (T, E) 
analysis. Using these definitions [Komatsu. Spergel & Wande ill d2005h define the one-point mixed-skewness involving the fields A 1 (r, Cl) and B J (r, Cl) 

(i,j,keT,E): 
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rdr / dflA u \r,Q)B v (r,Q.)B w (r, fi). (50) 



sA u B v B w _ 

S A B B can be used to estimate f l £ L , but such a radical compression of the data into a single number restricts the ability to estimate contamination 
of the estimator by other sources of non-Gaussianity. As a consequence, we construct a less radical compression, to a function of I, which can be used 
to estimate /]v£, but which can also be analyzed for contamination by, for example, foregrounds. We do this by constructing the integrated cross-power 
spectrum of the maps A u (r, Cl) and B v (r, Cl)B w (r, Cl). Expanding B 2 in spherical harmonics gives 

[B v (Cl,r)(r)B w (Cl,r)] lm = J d(lB v (r, ti)B w (r, Cl)Y lm (Cl) = £ £ (r) ( ^ ^ ) Iun"af, m ,a^ m „ (51) 

I'm' I" m" 

and we define the cross-power spectrum C Z A,S (r) at a radial distance r as 



C ^ v B w (r) = _l_^Real{AL(r)[B v (r)S w (r)] !m } I 

m 

Integrating over r we find: 



(52) 



(53) 



This integrated cross-power-spectrum of B 2 (r, tl) and A(r, O) carries information about the underlying bispectrum B wl /i, as follows: 
" b v b w 1 V^V^ , / / «' 



°' 21 + 



TlEE E J "'<" ( I L> l l» ) J r 2 dr{^(r)^(r)A^(r)}aLa^a^ m „. (54) 

m i'm' I" m-" ^ ^ 



C{ 21 + 



Similarly we can construct the cross-power-spectrum of the product map A B (r, Cl) and B (r, f2), which we denote as C, 

» fl v. fl w ^ /„„„ ( ^ ^ ^ ) / r 2 dr{A w (r)^(r)^ W }aLa^a- m „. (55) 

m Cm' l"m" 

Using these expressions, and the following relation, we can write this more compactly in terms of the estimated CMB bispectrum 

lVl " ~ I m m ' m " I a lmO-l'm'0-l"m" 

mm' m" 

from which we compute our new statistic, the bispectmm-related power spectrum, C\ oc as 

^" B v BW ^ ( ^, b v bW + . aMb v, b w + ^ fl w, BV) = £ S^f^lcTlC-riC- 1 ]^^ •} 07) 

W'V'W !' i" 

T TT EE S E^^ 

Clearly, in a joint analysis, the /rare skew-spectrum such as C ; ' or C l ' discussed in the previous section gets generalized to C, ' , etc, and 

the construction of inverse covariance weighted fields mixes temperature and i5-type polarisation. Hence, each of these estimators carries information 
from all possible types of mixed bispectra. If we combine various estimates into a unique power spectrum it compresses all available information for T 
and _E-type polarisation maps: 

6 ( 2 ,1) = cfB-B^ (5g) 

where B\°f v , is the bispectrum for the local }nl model, normalized t o i £^£ = 1. We can now use s t andard statistical techniques to estimate Note 
that if we sum over all / values then we recover the estimator S pr im of lKomatsu. Spergel & Wandelj J2005I) . which is the cross-skewness of ABB: 

S^ U V W =2^(21 + l)(Cf ' B B + Cf B ' B + Cf B ' B ). (59) 

i 

If we sum over all possible triplets ijk, we can recover 53 typically used in the literature and considered previously dMunshi & Heavensll2010l) . 

S 3 = S^ uBvBw . (60) 
uvw 

Though S3 compresses all available information at the level of the bispectrum in temperature and polarization maps, it clearly also has less power to 
distinguish any effect of systematics. These can be studied in more detail if we carry out individual estimates which break up the total into the estimates 
resulting from their various linear combinations. The method develops here a simple extension of previously used estimators both for one-point estimators 
as well as two-point estimators or the associated power spectra. Such methods will be useful with future surveys with higher signal to noise for polarization 

measurements. 

Recently, Calabrese et al. (2009) studied the impact of secondaries in estimation of primordial non-Gaussianity using temperature data. Further 
studies by Hika ge et alj j2009h used Fisher analysis to study joint two-to-one analysis. The estimator introduced above shows how individual linear 
combination of mixed bispectra can also be used for joint estimation especially while dealing with non-primordial contamination. 



12 D. Munshi et al. 



4.2 Power Spectra related to the Trispectrum 

In a recent paper, iMunshi et ai] J2009I) extended the earlier studies by Munshi & Heavens (2010) to the power spectrum related to the trispectrum in 
an optimal way. This statistic was applied to WMAP 5-year data in (Smi dt et al.ll2010h to provide the first constraints on tnl and jnl, the third order 
corrections to primordial perturbations in non-Gaussian models. The PCL estimators that we considered before are generalizations of these estimators 
from the temperature-only case to a fully joint temperature and polarization analysis. For a given choice of mixed bispectrum, a pair of corresponding 
power spectra can be defined which we optimize in this section for arbitrary partial sky coverage and instrumental noise. 



(3 1) 

4.2.1 Estimator for K.\ ' 

(31) (2 2) 

Moving beyond the bispectrum, we can construct estimators of the two power spectra we discussed above, C\ ' and C) ' . We show how to decompose 
these entire estimators to various choice of mixed bispectrum and compress the information optimally to define an unique estimator for each power spectra. 
The optimized versions for /C, can be constructed by cross-correlating the fields A u (ri, tl)B v (ri, Cl)B w (r2, ft) with B x (r-2, ft). In the first case 
the harmonics depend on two radial distances (n, r 2 ) for any given angular direction. For a specific combination of U, V, W and X which we can choose 
either to be temperature T or i?-type polarization E we can define a corresponding estimators. We will eventually combine all various contributions that 
we recover from these combinations to define a single estimator /C, 3 ' 1 ' which will generalize the estimator introduced in iMunshi & H eavens 201(3) for 
analysis of temperature-only data. 

A u (r 1 )B v (r 1 )\ tm = / A u ( ri ,Cl)B v fa, ft) Y^ft) dft; B(r 2 ) X \i m = J B x (r 2 , ft) F^(ft) dft. (61) 

Next, we construct the field C uv (ri, r 2 ) = Xw m -^H r i> r 2)A u (r{)B v (ri)\i m Yi m . If we now form the product of C uv (r±, r 2 ) and A w (r2) and denote 
it as D uvw (ri, r 2 ) = C uv {r\, r 2 )A (r 2 ) the spherical harmonic transform of this product field is represented as DY^ W (n, r 2 ). Finally we compute 
the cross-power spectra between D wvvv (ri,r2) andB x (r2). Wedenoteitby J t A B A - B (n, r 2 ) which also depend on both radial distances n and 

r 2 : 

jAU S v A",B* (rijra) = _l_ £ Re£j [ { ^ {ri ,r 2 )} lm {B x (r 2 )}l m ] . (62) 

m 

The construction for the second term is very similar. We start by decomposing the real space product B u (r, Cl)B v (r, Cl)B w (r, ft) and M (r, ft) in 
harmonic space. There is only one radial distance involved in both of these terms. 

B u (r,Cl)B v (r,Cl)B w (r,Cl)\ lm . = / [B u (r,ft)B v (r,Cl)B w '(r,ft)] F ; * m (ft) dft; M x (r,Cl)\ lm = [ M X (r, ft^^ft) dft. (63) 



Finally, the line-of-sight integral which involves two overlapping contributions through the weighting kernels for the first term and only one for the second 
gives us the required estimator: 



K yvw,X) =4/ 2 / r 2 dri J r l dr2 j«- «" "■•»~(r u r a ) + 2g, a J r'drC"" H ' a " ^ (r). (64) 

Next we show that the construction described above does reduces to an optimum estimator for the power spectrum associated with the trispectrum. 
The harmonics associated with the product field B u (n) B v (n) B w (r 2 ) can be expressed in terms of various f3(r) functions: 

B [ri)B (n)B (r 2 )| !m = a iimi a i 2 m2 a ! 3 m 3 a h (ri^fa)^ M^iU G Lh i ■ (65) 

The cross-power spectra J r ; ASA,S (ri, r 2 ) can be simplified in terms of the following expression: 

Jl (n, r 2 ) = ^— - 2^( _1 ) | i? i( r i. ? '2)a il (ri)/3 i2 (ri)a i2 (r 2 )/J; (r 2 )j {a hrni a hrn2 a l3m3 a lrn )g h l 2L 2 G h f L . (66) 

LA/ m 

The Gaunt integral describing the integral involving three spherical harmonics is defined as follows: 



rmi m 2 m 3 _ / (2il + l)(2k + l)(2i 3 + l) h h h \ h h l 3 \ (f ~ 

y hhh - V 4?r ^ o J \ mi m 2 m 3 ^ ' 1 J 

The second terms can be treated in an analogous way and the result takes the following form: 

r B u B V M W ,B X I \ 1 \ ^ i\A/V^ f oU I \qV ( \ W/ \qX{ vl l~U ~V ~W ~X \ n m. x m 2 M n m 3 mM ,an 

A ( V > = 2lTll^( > Z {AitO^^lOf^ WA W) (ai 1 m 1 ai 2ra2 ffli 3m3 ai m )6 lll 1 2i 2 G t J L ■ (68) 

LM m 

Finally, when combined these terms as in Eg. J64b. we recover the following expression: 



p.(UVW,X) _ 1 \ ^ \ " 1 rp-liWU'rp-liVV' rp-liWW'rp-li^'^liV^ rri 

; - 2F+T 2^ 2^ (2L + 1) 1 Jil [ Ji2 1 J ' 3 [ Ji w < 3 ^ [Ll 



i J / j /;:. /. 



(69) 
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We have subtracted the Gaussian component from the estimator in the last step. This is done by simulating Gaussian maps in a Monte-Carlo chain and 
using the same mask and the noise as the real data, weighting functions used for real data are also used on the Gaussian realisations. The ensemble average 
of the Gaussian realisations are then subtracted from the estimates from the real data. 

We can sum over all possible mixed bispectra to construct the following combined estimator: 

^ 3,1) = E £l" VW > X) . (70) 
uvwx 

The estimator ]C\ UVW > X > depends linearly both on f^ L and <?nl. In principle, we can use the estimate of /nl from a bispectrum analysis as a prior, or 
we can use the est imators 5; , /Cj and /C; 3 ' 1 ' to put joint constraints on /jvl and gNL- The former is better from a signal-to-noise point of view 
JSmidt et alj201Ch . Computational evaluation of either of the power spectra clearly will be more involved as a double integral corresponding to two radial 
directions needs to be evaluated. Given the low signal-to-noise associated with these power spectra, binning will be essential. 

4.2.2 Estimator for ^ uv ' wx '> 

In an analogous way the other power-spectra associated with the trispectrum can be optimized by the following construction. We start by taking the 
harmonic transform of the product field A(r, Q)B(r, Q) evaluated at the same line-of-sight distance r: 

A u (r,Cl)B v (r,U)\ lm = / A u (r)B v (r) Y t * m (n) dCl; B w (r,Cl)M x (r,Cl)\im = / B w {r)M x {r) Y^il) dtl, (71) 



and contract it with its counterpart at a different distance. The corresponding power spectrum (which is a function of these two line-of-sight distances n 
and r 2 ) has a first term 

jA" b^.a™ B*^^ = J r I ( riir2 )A w (ri ) n)B v (ri ) n)| Im .A(r 2) n) w B(r2,n) ;f |r m (72) 

m 

Similarly, the second part of the contribution can be constructed by cross-correlating the product of 3D fields B (Cl, ri)B v (Cl, n) against 
B w (ft, r 2 )M x (Q, r 2 ) evaluated at the same radial distance r 

cf « fl v, B w M * (r) = _1_ £ bU ^ A)sV(r) h )\ lm B w {r, h)M x {r, n)|f m . (73) 

m 

Finally, the estimator is constructed by integrating along the line of sight distances: 

&y V,W * ) =*&J&r 1 J ridr 2 J l AUBV > AWB \r ll r 2 ) + 2 gnl J r*dr Cf M * v .* w "* (r). (74) 
To see they do correspond to an optimum estimator we use the harmonic expansions and follow the same procedure outlined before: 



r B u B V ,B W M X 



/ \ 1 \ ^ / i\m \ ^ f qU { \ qV r \ qW / \ X / \°\ / ~U ~V ~W ~ X \ ^m-imom ^m^mAm 



(75) 



rrA u B v ,A W B x I \ 1 V~V iW V f ri/ \ Ui \ qV i \ Wi \ qX / \1 i -U -V ~W -X i/imimim^tnsmim m 

Jl (ri,r 2 ) = 2i^lZ^(- X ) \ F '('-i ! r2)Q il (ri)/3 !2 (r 1 )a i3 (r 2 )/3, 4 (r 2 ) j {a hmi a hm2 a l3m3 ai im4 }g h l 2l 2 Q h f it -(76) 

Here we notice that j^b,ab ^ j s mvar j ant unc jer exchange of ri and r 2 but j BB - BM ( ri r2 ) j s no t Finally, joining the various contributions to 
construct the final estimator, as given in Eq.J74t. which involves a line-of-sight integration: 

U'V'W'X' ijmj 34 
and the combined estimator 

jC (a, 3 ) = ^( W v,w*)_ (7g) 

uvwx 

f 2 2 ) (3 1 ) 

The prefactors associated with / NL and jnl are different in the linear combinations K.\ ' and K.) ' , and hence even without using information from 
third-order we can estimate both from fourth order alone. 

Similarly, the one-point cumulant involving both temperature and i5-type polarization at fourth order can be written in terms of the the mixed trispectra 
as follows: 

= J] (2Z + l)£p< 2) =5^(21 + 1) Y, ^ (WV ' VV "°=E( 2Z + 1 )^ (3 ' 1) =E( 2Z + 1 ) E ^ VW,X) - (79) 
1 1 uvwx 1 1 uvwx 

It is also possible to carry out the sum over the harmonics I without summing over the field types. In this case we recover independent one point estimators 
associated with each type of mixed trispectrum. 
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K uvwx = 5^(21 + i)t^ v ' wx) = + l)tj uvw ' x \ (80) 
i i 

We will next consider the correction terms for these estimators due to absence of absence of spherical symmetry - which may be broken either because of 
mask or due to the presence of detector noise in an arbitrary scanning strategy. 

4.3 Correction in the Absence of Spherical Symmetry 

It was pointed out in Babich ( 20051) : ICreminelh et al.l {2006); Yadav et al. ( 2008) that in the presence of a partial sky coverage, e.g. due to the presence of 
a mask or because of the galactic foregrounds and the bright point sources, as well as, in the case of non-uniform noise, spherical symmetry is destroyed. 
The estimators introduced above will then have to be modified by adding terms which are linear in the observed map. The corrective terms are incorporated 
using Monte-Carlo techniques. A more general treatment which involves computationally expensive inverse covariance weighting will be discussed later. 
The treatment discussed here is nearly optimal though ignores mode-mode coupling dominant at low I, 



4.3.1 Corrective terms for C, ' 

The (linear) corrective terms are constructed from correlating the Monte-Carlo (MC) averag ed (A u (r,Cl) B v (r,Cl)) alm product maps with the input 
B w (r, Cl) map. The mask and the noise that are used in constructing the Monte-Carlo averaged product map are exactly same as the observed maps and 
the ones derived from them such as A or B. Mode-mode coupling is important at low angular modes, and we consider the full case later, but for higher 
frequency modes, we can approximate the linear correction to the local shape: 
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(81) 



where f s k y is the observed sky fraction. 

The C;s such as C( A describe the cross-power spectra associated with Monte-Carlo (MC) averaged product maps {A(n, r)B(n, r)} constructed 
with the same mask and the noise model as the the observed map B. Likewise, the term C^ B ' B ' > denotes the average cross-correlation computed from 
MC averaging, of t he product map construc ted from the observed map r) multiplied with a MC realization of map B(Q, r) against the same MC 
realization B(0., r). ICreminelli et aD {2006) showed via numerical analysis that the linear terms are less important in the equilateral case than in the local 
model. The use of such Monte-Carlo maps to model the effect of mask and noise greatly improves the speed compared to full bispectrum analysis. 

The use of l i near t erms was found to greatly reduce the scatter of this estimator, thereby improving its optimality. The estimator was used in 
lYadav & Wandeltj d2008h also to compute the }nl from combined T and E maps. The analysis presented above is approximate, because it uses a crude 
f s ky approximation to deconvolve the estimated power spectrum to compare with analytical prediction. A more accurate analysis should take into account 
the mode-mode coupling which can dominate at low I. The speed of this analysis depends on how fast we can ge nerate non-Gaussian map s. The general 
expression which includes the mode-mode coupling will be presented in the next section. However it was found bv lYadav & Wandeltl ( 120081) that removing 
low Is from the analysis can be efficient way to bypass the mode-mode coupling. A complete numerical treatment for the case of two-point statistics such 
as C, ' B will be presented elsewhere. 
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4.3.2 Corrective terms for the Estimator fC t 
For the four-point terms, we need to subtract linear and quadratic terms: 
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The expressions are similar for the other terms that depend on only one radial distance 
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To simplify the presentation we have used the symbol A u (ri, Cl) = A^ 1 ; A u (r, Cl) = A u and so on. Essentially we can see that there are terms which 
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are linear in the input harmonics and terms which are quadratic in the input harmonics. The terms which are linear are also proportional to the bispectrum 
of the remaining 3D fields which are being averaged. On the other hand the prefactors for quadratic terms are 3D correlation functions of the remaining 
two fields. Finally putting all of these expressions we can write: 



C.(UVW,X) .,2 I 2, 1 2, j-A u B v A w ,B x , \ , r, I 2, ~B U B v M vv ,B" 

K-i — 4/nl / r 1 dn / r 2 dr 2 J l {n, r 2 ) + 2g NL I r dr£ l 



(88) 

From a computational point of view clearly the overlap integral Fi,{ri,r2) will be expensive and may determine to what resolution ultimately these 
direct techniques can be implemented. Use of these techniques directly involving Monte-Carlo numerical simulations will be dealt with in a separate 
paper (Smidt et al. in preparation). To what extent the linear and quadratic terms are important in each of these contributions can only be decided by 
testing against simulation. 



4.3.3 Corrective terms for the Estimator K. l 

The unbiased estimator for the other estimator can be constructed in a similar manner. As before there are terms which are quadratic in input harmonics 
with a prefactor proportional to terms involving cross-correlation or variance of various combinations of 3D fields and there will be linear terms (linear 
in input harmonics) with a prefactor proportional to bispectrum associated with various 3D fields. 
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The terms such as K, l ' (n 5 T2) can be constructed in a very similar way. We display the term /C z ' (ri , r2) with all its corrections included. 
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The importance of the linear terms depends greatly on the target model being considered. For example, while linear terms for bispectral analysis can 
greatly reduce the amount of scatter in the estimator for local non-Gaussianity, the linear term is less important in modeling the equilateral model. In 
any case the use of such Monte-Carlo (MC) maps is known to reduce the scatter and can greatly simplify the estimation of non-Gaussianity. This can be 
useful, as fully optimal analysis with inverse-variance weighting, which treats mode-mode coupling completely, may only be possible on low-resolution 
maps. 



p.(UV,WX) At 2 12, I , ^A" B v ,A VV B x , , . „ I 2, ;B U B v ,B yv M rt , \ 

M —4/nl / r x dr\ I r2dr2j l (ri, r2) + 2qnl I r dr£ [ (r). 



(95) 



The corrections to one-point estimators can be recovered by performing appropriate sums. In the next section we use direct summations and proper 
modelling of the covariance matrix as opposed to the Monte-Carlo techniques used here. However in certain situation it may be difficult to model the 
covariance matrix in an accurate way, we also provide analytical results which can handle such situations. 



5 EXACT ANALYSIS OF OPTIMAL ESTIMATORS 

In the previous section we relied on Monte-Carlo simulations to model the effect of finite sky coverage, noise as well as other observational artefacts. 
The resulting method nearly optimal and for most cases where mode mode coupling is not to strong i.e. for near all-sky coverage can be efficient. To take 
mode-mode coupling properly into account which is the case for low multipoles we need to model the covariance matrix as accurately as we can. In this 
section we tackle the case where inverse covariance weighting is required. The resulting method is optimal and can provide accurate results for studies 
with degraded resolution maps where cross-contamination to primordial non-Gaussianities coming from secondaries is minimum. 



5.1 The Power Spectrum Related to the Bispectrum 

The general expression for the bispectrum estimator was developed bv lBabicf] d2005h for arbitrary sky coverage and inhomogeneous noise. The estimator 
includes a cubic term, which by matched-filtering maximizes the response for a specific type of input map bispectrum. The speed of this analysis depends 
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on how fast we can generate non-Gaussian maps. The linear terms vanish in the absence of anisotropy but should be included for realistic noise to reduce 
the scatter in the estimates; see Babich (2005) for details. We define the optimal estimator as: 



E X L yZ [a] 



=E^w[g E E E *2B?(£ I i>) 

" X'y'Z' MM' ll'limm'm. 



~ l^lm,l'm'] ^ '(Pi'M'.ljmjl a i 2 ™2 ) ~~ Pi. Af.im] ( Pi'm' ,l 2 m 2 a h™2 ) _ pLM,!m]' > ' (Pl'm'^tna ] a i 2 m 2 )} ! 

A", y, 2, x',y',z' g {t, s} (96) 

where iV^/ is a normalization to be discussed later. A factor of 1/(21 + 1) can be introduced with the sum ^2 M , if we choose not to introduce 
the N LL i normalization constant. This will make the estimator equivalent to the one introduced in the previous section. As the data is weighted by 
C _1 — (S + N)^ 1 , or the inverse covariance matrix, the speed of this analysis depends on how quickly we can generate non-Gaussian maps, addition 
of higher modes will reduce the variance of the estimator. In contrast, the performance of sub-optimal estimators can degrade with resolution, due to the 
presence of inhomogeneous noise or a galactic mask. However, an incorrect noise covariance matrix can not only make the estimator sub-optimal but it 
will make the estimator biased too. The noise model will depend on the specific survey scan strategy. Numerical implementation of such inverse-variance 
weigh ting or multiplication of a map by C _1 can be carried out by conjugate gradient inversion techniques. Taking clues from Smith & Zaldarriaga 
( 2006), we extend their estimators for the case of the skew spectrum. We will be closely following their notation whenever possible. First, we define 
Q l [a] and its derivative dimQh [a] . The required input harmonics ai m are denoted as a. 

M l'm',l"m" ^ ' 

*x /\x.yz, ,_ 1, v~~* r,xyz I L I' I' \ y z , no * 
dimQ L [a] = -S L i > B un „ \ , „ o,/ m ,oi// TO «; (98) 



V m' .1" m" 



m m m 



M I'm' * ' 

These expressions differ from those for the one-point estimators by the absence of an extra summation index. If summed over the free index L, the 
expression for Ql reduces to a one-point estimator. Q x ' yz [a] represents a map as well as di m Q X [a], however Q x ' yz [a] is cubic in input maps 
a; m s where as the derivatives di m Q x ' y2 [a] are quadratic in input. The expression for the derivative, is different, when the derivative is taken w.r.t the 
field (e.g. X in this case) associated with the free indices than when it is taken with respect to a field y or Z whose indices are summed over. 
The skew-spectrum can then be written as (the summation convention is assumed for the next two equations): 



(100) 



i%' y *[a] = [n- 1 ]^ ^Q x L ; yz [C' 1 a]^^[c- 1 a\L(dLQ x j yz [c-\^)Mc)^ ; s e (x,y,z) 

Here ()mc denotes the Monte-Carlo averages. The inverse covariance matrix in the harmonic domain [C ,xv ] ; ~;J ni i 2Trl2 = { a i" 1 m 1 a f m } _1 encodes the 
effects of noise and the mask. For all-sky and in the signal-only limit, it reduces to the usual [C ,_1 ]*^ i l2Tn2 = xy 5u'&mm' The normalization of the 
estimator which ensures unit response can be written as: 

N LV =lYl [<{^ mi QL[C" 1 a]}[C' Iimill2m2 ] SS '{C ! Q^[C- 1 fl]}) 
SS' 

-{<^ mi 0i[C- 1 a])}K-i, lil!ima ] ss '{<< mB Q £ ,[C- 1 o])}]; Q L =Q x ' yz ; S,S'eX,y,Z. (101) 

In the above expression, and in those that follow, we will not explicitly display the superscript on the normalization matrix N LL i and Ql as they are 
obvious from the context. Summing over repeated indices is assumed, and the second term ensures subtraction of terms with self-coupling(s). We will be 
using the following identity in our derivation: 

^Si,l 2 m2 = V^Jimi^iT^) ~ (P a ]l 1 m 1 [C <^ 2 m 2 ) = ^""^ \P \hmx ,l a *n a ^Um a , l b m b [C ]hm 2 ,lb^b' (102) 

X'y' l a ™>a lb m b 

The Fisher matrix, encapsulating the errors and covariances on the El, for a general survey associated with a specific form of bispectrum can finally be 
written as: 

F W = ^ ( (1) a££, + (2) aff, + w affi + (3) a?« + (s) ag« + (4) ag«) . (103) 

Using the following expressions, which are extensions of lSmith & Zaldarriaga (2006), we find that the Fisher matrix can be written as a sum of two a 
terms a pp and eft®. The terms involved a correspond to coupling only of modes that appear in different 3j symbols. Self-couplings are represented 
by the beta terms. The subscripts describes the coupling of various I and L indices. The subscript PP correspond to coupling of free indices,i.e one free 
index L\ with another free index L2 and similar coupling for indices that are summed over such as h, I2 etc. Similarly for subscript QQ the free indices 
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are coupled with summed indices. Couplings are represented by the inverse covariance matrices in the harmonic domain e.g. C, * LM denotes coupling 
of mode LM with Im. 

„pp _ n xyz n xyz ( Li h l[ \ I L 2 h 1' 2 \ r /r ] xx tr ^y\r \ zz 

W a L 1 L 2 - 2^ 2^ B L lh l^L 2h l' 2 I Ml mi m / ) [ M2 m2 m < 2 I [U lMu L 2 M 2 \ lC !imilIama j IC,/^,,^] 
M!,M 2 lil'.mim'. V 7 V 7 

_QQ _ R .WZ r j;jz f h 'i \ I L 2 h 1' 2 \ ^ \ x y\r ^ x \r \ zz 

W a L lL2 - 2^ 2^ B l 1 i 1 v 1 b l 2 i 2 v 2 I Mj mi mi M M a m 2 m 2 J L c ^i^i^ 2 -2J N iroi ,£ 2 M 2 j L^ S ,i< S J 

a£f i2 = (■ LxjipCLajap ; o£« a = ( Lijip(.Lajaj a) . (104) 

Similar results for ( 2 ) <^l^l 2 (2-4) a 2^L 2 can be obtained from permutative reordering of the above results. The ordering of the multipole indices and 
that of corresponding fields denoted by X is important. For each different choice of field triplets we will have a different set of skew-spectra associated 
with the bispectrum. If we choose to have the same triplets for the primed and unprimed fields then we will recover the Fisher matrix associated with that 
specific choice. However, if we decide to choose a different set of triplets then we will have the information about the level of cross-contamination from 
one type of power spectra to another. Within a specific choice of triplet e.g. TEE choice to associate the free index L to a given field type e.g. T or E 
will generate two different skew spectra from the same bispectrum B TEE . Results presented above will reduce to those of ISmith & Zaldarriaga (2006) 
when further summations over free indices L\ and L 2 are introduced to collapse the two-point object to the corresponding one-point quantity. The /3 
terms that denote cr oss-coupling are not presen ted here as they do not appear in the final expressions for the Fisher matrix. A detailed analysis of these 
terms is presented in lMunshi & Heavens! ( 120 101) and can be extended in a very similar manner. If we sum over LL' the Fisher matr ix reduces to a scalar 
F = V] , F LL i with, = a-12' ~ a anc ^ Pll' = Pll' = Pll' ~ P' wrlere Q ' P an( ^ F are exactly the same as introduced in lSmith & Zaldarriagal 
(2006)" for one-point estimators. 



5.1.1 Joint Estimation of multiple Bispectrum-related Power-Spectra 

The estimation technique described above can be generalized to cover the bispectrum-related power spectrum associated with different sets of bispectra 
(X,Y), where X and can Y can be one of the combinations from the set { (TTT), (TTE), (TEE), (EEE)} Loc/eq . In addition to considering triplets 
corresponding to a given primordial bispectrum, we can as well consider joint estimation of different initial primordial bispectrum such as local (loc) 
bispectrum or equilateral (eq). 

E r L [a] = ^[F- 1 ]^'' j^'P^l ~ X)[C^ 1 a]to.(^Q^[C _1 a]>«c)| ; T, T' G (TTT), (TTE), (TEE), (EEE) loc/eq ; SeT,E. (105) 

The associated Fisher matrix now will consist of sectors F^,,FY\i and F^S ■ The sector TV and T'T' will in general will be related to errors associated 
with estimation of bispectra of F and F' types, whereas the sector rr' will correspond to their cross-correlation. Clearly with a given mixed bispectra 
type it is possible to have different estimators by associating a specific type of field T or E with the index which is not summed over. 



(106) 



Here the dots represent contribution from other terms represented by iCt, which are obtained from permutation of various indices. We have introduced the 
following notation above: 

iK£<] rr =J2 Bl hl{ Bl, hl , 2 ( ^ ^ ^/)( M / ^ 1^\[Clm,l>m'] XX [Ci imi ,i 2 m 2 ] yy [C,/ mi ,,/ m /] M (107) 

MM' Ijl'.miTn'. V 7 V 7 

and a similar expression holds for the other [a££/] terms as well as the [a^v] terms. We can also use the above formalism to obtain the cross-correlation 
of a given skew-spectrum type from different kinds of primordial bispectrum, as well as, say, local type and equilateral types of non-Gaussianity. 



5.1.2 All-sky Homogeneous Noise 

The above expressions are very general results for arbitrary sky coverage due to a specific scanning strategy. Our approach can deal with complications 
resulting from partial sky coverage and inhomogeneous gaussian noise. Any residual non-gaussian noise will have to be subtracted out and will need 
more elaborate analysis of variance estimation. 

If we now take the limiting case when we have all sky coverage and homogeneous noise we can recover analytical results which are useful for 
comparing various planned and ongoing surveys. In the the all-sky limit, the covariance matrices are determined entirely by signal and noise power 
spectra. To simplify the general expressions derived so far for the case of all-sky coverage and uniform noise we will use the following expression: 

xy l ( {/d? T l/d TE \ 

Cl 1 m 1 ,l 2 m 2 = [C lil^l^mima — I \j^ E l/df E J ^ l l h 2^ m l m 2- (108) 



We recover following expression for the case of temperature: 
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If we assume that there is no correlation between the temperature and _E-type polarization then df = Cf which reduces to the temperature-only result. 
Notice that in this case the corresponding ;a pp functions and iof 1 ^ functions become degenerate. Expressions for the case of estimation error of 

T EE E TE 

E L ' L , , E L ' L , can be obtained using the following expressions for the related Fisher matrices: 
i f \n TEE ^ ln TEE ] 2 1 

F " B = M \ Sll ' D b "*J 2 (dfdfdf + Iffdj) + D B "™ (djdjdf + IfdJiJ + dppf + Sfdpf) } <lll) 

The Fisher matrices for F®}J T and F^f E can also be constructed in a similar manner. 

Using a specific form for the bispectrum bj°j 2i , such as the local model, the Fisher matrix elements can be further expressed in terms of the transfer 
functions a x and /3 Y and the associated power spectra C XY using the definitions of iOt pp and icfi® introduced before: 

^r- ^( 2i+ i„ 2£ ' +1 ,B a+ i)(S o oJwkw 

x |y r 2 dr(aZ(r)ft,(r)p?(r) + al(r)pUr)P?(r) + af(r)tf(r)Pl(r))^ ; T = Xyz (112) 
1 [^] IT = %(2L + l) E (2Z + l)( 2 r + l)( J J J ) 2 -^ 

ii' v / L l V 

xy\- 2 dr(af(r)^(r)l3^(r) + af(r)^(r)^(r) + a^(r)l3f(r)^(r)^ ; r = tfyZ. (113) 

The power spectra Ci appearing in the denominator take contributions from both the pure signal (i.e. the CMB) and the detector noise. It is possible to bin 
the estimates in sufficiently large bins that these are practically uncorrelated estimates for an experiment such as Planck with very high sky-coverage. A 
joint analysis will combine the the results from all possible estimators. The cross-correlation among various estimators (characterized by different choice 
of Xy and Z) can be computed following the same techniques. A detailed analysis of the singularity structure of the error-covariance matrix will be 
presented elsewhere. 

5.2 Power Spectra related to Trispectra 

Extending our analysis to the case of four-point correlation functions involving both temperature and polarization data, we will next consider the case of 
one-point and two-point estimators which are related to the mixed trispectra. 

5. 2. 1 One-point Estimators 

We will use an inverse- variance weighting for harmonics recovered from the sky. The covariance matrix, expressed in the harmonic domain, {aY m a^, m , ) — 
[C _1 ]^ ( , m , is used to filter out modes recovered directly from the sky a; m . We use these harmonics to construct optimal estimators. For all-sky coverage 
and homogeneous noise, we can we recover the results derived in the previous section. We start by keeping in mind that the trispectrum can be expressed 
in terms of the harmonic transforms a; m , which can either be temperature multipoles or polarization multipoles : 

4V( 2! + 1 'EEH) M (i t b t _ L M )at ma ...af dmd (iea, b ,c, d ), (U,V,W, X) e T, E. (114) 

m. t M V 7 V 7 

Based on this expression we can devise a one-point estimator. In the following discussion, the relevant harmonics can be based on partial sky coverage. 
Q «vw. [a] ly H) My A(i i)r yJ 1« h L \f l c h L \ a u ... a ? (115) 

LM li-rrii v 7 v 7 

The term A(k, L) is introduced here to avoid contributions from Gaussian or disconnected contributions: A(h, L) vanishes if any pair of ks becomes 
equal or L — which effectively reduces the trispectra to a product of two power spectra (i.e. disconnected Gaussian pieces). Its value is unity for the 
connected terms. We will also need the first-order and second-order derivative with respect to the input harmonics. The linear terms are proportional to 
the first derivatives, and the quadratic terms are proportional to second derivatives, of the function Q[a], which is quartic in the input harmonics. 

d lm Q W = 4t2^(- 1 ) l^ A ( l '' L ) T m c x ld [ m ma M )[m b m c -M J °«.»*<W (116) 
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^q( uvwx ) j a ] i s cubic in the input maps and the second-order derivative is quadratic in input maps (in terms of 
^(uvwx) ^ w hich is simply a number, these objects represent maps constructed from harmonics of 
the observed maps. The quadratic terms contribute only to disconnected parts and hence will not be considered. The optimal estimator for the one-point 
cumulant can therefore be written as follows. This is optimal in the presence of partial sky coverage and most general inhomogeneous noise: 

E uvwx [a] = A !q« vwx [C-\] -^[C-^UdLQr™*^])} ; S e (U,V,W,X). (117) 

I S,lm ) 

The terms which are subtracted out can be linear or quadratic in input harmonics. The linear term is similar to the one which is used for bispectrum 
estimation, whereas the quadratic terms correspond to the disconnected contributions and will vanish identically as we have designed our estimators in 
such a way that it will not take any contribution from disconnected Gaussian terms. The Fisher matrix reduces to a number which we have used for 
normalization. The direct summation over various harmonics as described above can be expensive computationally and determines to what resolution the 
numerical calculation can be performed. The input [C~ a] denote the entire set of inverse variance weighted harmonics for the Q and dQ. It is interesting 
to note that modes after inverse variance weighting are no longer pure and are linear combinations of Temperature and i3-type polarization modes. This 
true for the case of estimators too. A given estimator p uvwx though correspond to a specific choice of trispectrum takes contributions from modes 
which are themselves linear combinations of pure mode types. 

^ VW " = ^=(l) 2 E EE E E (-1) M (-1) M 'A(^;L)A(^L')^(L)7^^(L) 

U'V'W'X 1 LAI h'W (all Im) (all I'm') <= d 

h L \ f la Id L \( l' a l' b L' \ f l' c l' d L' 
m b M ) \ m c m d -M ) I m' m' b M' I I m' c m' d —M' 



l a m a ,l' a m'J •••Rjm j ,l>d +cyc.permf. (118) 



The cyclic permutations in these terms will include covariances involving all-possible permutations of the four fields involved in construction of the mixed 
trispectrum and the related power spectra. The ensemble a verage of this one-point est i mator will be a linear com bination of parameters / NL and <7nl. 
Estimators constructed at the level of three-point cumulants dSmith & Zaldarriagall2006l : lMunshi & Heavensll2010h can be used jointly with this estimator 
to put independent constraints separately on /nl and jnl- As discussed before, while the one-point estimator has the advantage of higher signal-to-noise, 
such estimators are not immune to contributions from an unknown component which may not have cosmological origin, such as inadequate foreground 
separation. The study of these power spectra associated with bispectra or trispectra can be useful in this direction. Note that these direct estimators are 
computationally expensive due to the inversion and multiplication of large matrices, but can be implemented in low-resolution studies where primordial 
signals may be less contaminated by foreground contributions or secondaries. 

Combining all possible choices of mixed trispectra it is possible to introduce one single number which represents the entire information content 
regarding the trispectrum from Temperature and iJ-type polarization in the absence of B-modes. The corresponding estimators and the Fisher matrix 
takes the following form: 

E W _ e uvwx ; F (4) = V" pUVWX 

uvwx uvwx 

The choice of estimator is related to the level of compromise one is willing to made to increase the signal-to-noise at the expense of losing the crucial 
power to distinguish contributions from various contributions. Clearly it is also possible to design estimators by fixing a subset of all available indices 
representing the choice of T or E. 



5.2.2 Two-point Estimators 

Generalizing the above expressions for the case of the power spectrum associated with trispectrum, we introduce two power-spectra which we have 
discussed in previous section in the context of construction of nearly-optimal estimators. The information content in these power spectra are optimal, and 
when summed for over L we can recover the results of one-point estimators. 

MUV,WX)r i 1 V-*/ -,\M ST^ r r u la V l b f T \ (la lb L \ ( Ic Id L \ u X , nm 

Ql H = 1 ^(-l) l^ T mA {L) [m a m b M { m c m d _ M J ^ • • • <W (120) 

The derivatives at first order and second order are as series of maps (for each L) constructed from the harmonics of the observed sky. These are used in 
the construction of line ar and quadratic term s. We have retained the overall normalization factor so that our estimator reduces to the temperature-only 
estimator introduced in Muns hietal] ( l2009h . 

T iiTTli V 7 V 7 

We can construct the other estimator in a similar manner. To start with, we define the function q(U' vwx '> j a ] an( j construct its first and second derivatives. 
These are eventually used for construction of the estimator E^' vwx \a\. As we have seen, both of these estimators can be collapsed to a one-point 
estimator Q uvwx [ a J_ As before, the variable a here denotes input harmonics a l f m recovered from the noisy observed sky. 
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M ST i,™; V ' V 7 

The derivative terms will have two contributing terms corresponding to the derivative w.r.t. the free index {LAI} and the terms where indices are summed 
over e.g. {lm}, which is very similar to the results for the bispectrum analysis with the estimator Q * y ' z [a] . One major difference that needs to be taken 
into account is the subtraction of the Gaussian contribution. The function A(h, L) takes into account of this subtraction. 



r-U ^(U,VWX)r , r V~~* V~~* A /; r rr« t^-L V ; 2 lrr \ (I h S \ I lc Id S \ V W 
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Using these derivatives we can construct the estimators e^,vwx) an( j ^j(p' v > w '*0 . 



t-,(U,VWX) 
L 



NZb \Q^ VWX) [C-'a] - ^\C-\}L{dLQf, yWX \C-\])^ ; 3 G W, V, W, X (125) 

4 WV ' W "' = *Eb ^^[C-'a] -^[C-MLidLQT^P-'a]}] ; S eU,V,W,X. (126) 

where summation over L' is implied. The quadratic terms will vanish, as they contribute only to the disconnected part. The normalization constants are 
the Fisher matrix elements F LL i which can be expressed in terms of the target trispectrum T^i* (L) and inverse covariance matrices C _1 used for the 
construction of these estimators. The Fisher matrix for the estimator fi^' vwx ^ j e piP^vwx) ^ e expressed as: 



[n-\ l , = F %r x) = (±) 2 E E (-i) M (-i) A/ '[^(^[n5;5;(5')] 

ST,S'T' (all lm,Vm') c d 

wmw) x ( £ £ £ S T ) ( £ 4 J ) ( L M , * ) 

X I [Cl,M,L'Af'] WW lCl a m a ,l'a m 'a) ^ "H >'i m i ] l-^lamcl^m'J +••■ 

+[^M,!^m' £l ] WV [Ci a m a ,L'A/'] VW [C; b m i ,,;' b m'J WW [C; c m c ,i^ m ' c ] +... |. (127) 

The first set of terms can be recovered from the first term by permuting the multipole indices while still keeping the coupling of the free indices LL' intact. 
Similarly the second set of terms represented by . . . can be recovered from the second terms but considering only coupling between free indices and the 
one that are summed over. There will be a total of six terms of first type and eighteen of the second type. Similarly, for the other estimator E^f v ' wx ^ , the 
Fisher matrix F^^> wx > can be written as a function of the associated trispectrum and the covarian ce matrix of various modes . For further simplification 
of these expressions we need to make simplifying assumptions for a specific type of trispectra, see lMunshi & Heavens! d2010h for more details for such 
simplifications in the bispectrum. 



[*-w^ w ^) 2 E£ E E £ m)(m t -m) 

M M' (all lm) (all I'm') c d \ / V / 

X ( k < M> ) ( M' K -M' )A(W;L)A«;L') { [C lamM ] UU ' . .. fa^^]** + cyc.pern,} (128) 

Knowledge of the sky coverage and the noise characteristics resulting from a specific scanning strategy is needed for modelling of [C~ ff m , . We 
will discuss the impact of inaccurate modelling of the covariance matrix in the next section. The direct summation we have used for the construction of 
the Fisher matrix may not be feasible except for low resolution studies. However a hybrid method may be employed to combine the estimates from low- 
resolution maps using the exact method with estimates from higher resolution maps using other faster but sub-optimal techniques described in previous 
section. In certain situations when the data is noise-dominated further approximations can be made to simplify the implementation; a more detailed 
discussion will be presented elsewhere. It is possible to sum over all possible trispectra to recover the entire information content: 

#(1,3) _ y^ E (uywx)_ F (i,3) _ y^ F (u,vwx)_ ^(2,2) = y^ E (u,vwx)_ ^1,3) _ y^ F (u,vwx) (129) 

I / j I LL / j LL' I / j I LL / j LL 1 

uvwx uvwx uv,wx uvwx 

Summing over the free indices we recover the one-point estimators and the corresponding Fisher matrices: 
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Figure 4. Various Trispectrum related power-spectra are plotted as a function of angular scale. An ideal all-sky no-noise experimental set up was used for computing the 

power spectra. Three different trispectra were considered for the combinations TTTT, EEEE and TETE. The left panels correspond to the estimator (21 + lJ/Cj 3 ' 1 ' and 

(2 21 

the right panels correspond to (21 + 1)/C, ' . Upper panels correspond to gjvL = 1 ar| d the bottom panels correspond to /atl = !• Cosmological parameters correspond 



i that of WMAP7 analysis lLarson et alj2010h . See text for details. 
B W =^E^ = $>?' 2) ; j*) = £^> = £ F (M>. (130) 

( l LL' LL' 

So far we have assumed that the covariance matrix can be modelled accurately. In the next section we will discuss the impact of not knowing the covariance 
matrix accurately. We will show that though the estimators still remain unbiased but they no longer remain optimal. 

5.2.3 Approximation to exact C _1 weighting and non-optimal weighting 

If the covariance matrix is not accurately known which is most often the case due to the lack of exact beam or noise characteristics, as well as due to 
limitations on computer resources to model it to high accuracy, it can be approximated. An approximation R of C~ then acts as a regularization method. 
The corresponding generic estimator can then be expressed as: 

£f[a] = ^[F- 1 ]^|gf,[fla]-^[i?a]f m (9f m Qf,[^])|; Z e {(UV,WX), (U,VWX)}; SeU,V,W,X. (131) 

As before we have assumed sums over repeated indices and (■) denote Monte-Carlo (MC) averages. As is evident from the notation, the estimator above 
can be of type ft(u,vwx) Qr j^(uv,wx) p of ^ collapsed case can also be handled in a very similar manner. 

E z [a]=^El = ^[F- 1 ] LL , {Ql\Ra]-\Ra]L{dLQl[Ra])} Z £ {(WV, WX), (U, VWX)}; SeU,V,W,X. (132) 

L LL' 

We will drop the superscript Z for simplicity, but any conclusion drawn below will be valid for both specific cases i.e. Z £ {(WV, WX), (U, VWX)}. 
The normalization constant which acts also as inverse of associated Fisher matrix Flu can be written as: 

Ell' = ({El){E l ,)) - ((E L )){(E L ,)) = i ^ { (df m Q L [Ra] [Cr^fm'j'm'^m'Qi' \M) ~ (d? m Q ' i\Bn\)[C~ x ]f£ v w <$LQ l'[M)}\ 

S3' 

S, S' G (U, V, W, X). (133) 

The construction of E LL i is equivalent to the calculation presented for the case of R — C _1 . For the one-point estimator we similarly can write 
F R = r j , F^ L ,. The optimal weighting can be replaced by arbitrary weighting. As a special case we can also use no weighting at all R = I. This 
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reduces the cost of the estimator drastically. Although the estimator still remains unbiased the scatter however increases as the estimator is no longer 
optimal. Use of arbitrary weights makes the estimator equivalent to a PCL estimator discussed before. In certain circumstances the use of a fast method 
can be very useful before applying more robust and optimal techniques. 



5.2.4 Joint Estimation of Multiple Mixed Trispectra 

It may be of interest to estimate several trispectra jointly. The different sources of trispectra can be of all primordial type such as from "adiabatic" and 
"isothermal" perturbations. Such an estimation can explore the joint error budget on parameters involved from the same data-set. In such scenarios it is 
indeed important to construct a joint Fisher matrix which will take the form 

Bi[o] = y^y^[-F -1 ]i£''-Ei'[a]; T,r' € Adiabatic, Isothermal. (134) 

rr' ll> 

The estimator E^[a] is generic and it could be either ij' 3,1 ' or E ( - 2 ' 2 \ Here X and Y corresponds to different trispectra of type X and Y, these could 
be e.g. primordial trispectra from various inflationary scenarios. It is possible of course to do a joint estimation of primary and secondary trispectra. The 
off-diagonal blocks of the Fisher matrix will correspond to the cross-talk between various types of bispectra. Indeed, a principal component or generalized 
eigenmode analysis can be useful in finding how many independent components of such trispectra can be estimated from the data. 
The cross terms in the Fisher matrix elements will be of following type: 
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(135) 



The expression displayed above is valid only for E^ UV ' WX \ exactly similar results holds for the other estimator E^ u,vwx \ For X — Y we recover 
the results presented in previous section for independent estimates. As before we recover the usual result for one-point estimator for Q 4 from the Fisher 
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matrix of Q L ' or Q L ' , with corresponding estimator modified accordingly. 
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A joint estimation can provide clues to cross-contamination from different sources of trispectra. It also provide information about the level of 
degeneracy involved in such estimates. It is also possible to do joint estimation involving two different types of power spectra associated with trispectra 
or to include even the power spectrum associated with the bispectrum. The results presented in this section can be generalized to include such cases too. 



6 CONCLUSION 

The ongoing all-sky survey performed by the Planck satellite will complete mapping the CMB sky in unprecedented detail, covering a huge frequency 
range. It will provide high resolution temperature and polarization maps which will provide the cosmological community with the opportunity to constrain 
available theoretical models with unprecedented accuracy. Recent detection of non-Gaussianity in WMAP data has added momentum to non-Gaussianity 
studies. 

The temperature and polarization power spectra carry the bulk of the cosmological information, though many degenerate early universe scenarios can 
lead to similar power spectra. Higher-order multi-spectra can lift these degeneracies. The higher-order spectra are the harmonic transforms of multi-point 
correlation functions, which contain information that can be difficult to extract using conventional techniques. This is related to their complicated response 
to the inhomogeneous noise and partial sky coverage. Analysis of the higher order polarization statistics is more complicated not only by relatively low 
signal-to-noise ratio and their spinorial nature but also because of uncertainties in modeling the polarized foreground. A practical advance is to form 
collapsed two-point statistics, constructed from higher-order correlations, which can be extracted using conventional power- spectrum estimation methods. 
If the polarization data can be analyzed and handled properly it can help to further tighten the constraints on parameters describing non-gaussianity that 
are achieved by analys is of temperature dat a alone. 

In a recent study. iMunshi et al. (2009) studied and developed three different types of estimators which can be employed to analyze these power 
spectra associated with higher-order statistics such as bispectrum or trispectrum for analysis of temeprature data. In this work we include polarization on 
top of this to further tighten the constraints. The analysis of polarization data makes the analysis bit involved because of the spinorial nature of the data. 
We start with the MASTER-based approach dHivon et al.l2002l) which is typically employed to estimate pseudo-C;s from the masked sky in the presence 
of noise These are unbiased estimators but the associated variances and scatter can be estimated analytically with very few simplifying assumptions. 
We extend these estimators to study higher-order correlation functions associated mixed multispectra such as bispectrum and trispectrum which involve 
both temperature and polarization. Ignoring the contributions from B-polarization makes these estimators much simpler. These estimator can be very 
useful for testing simulation pipeline running Monte-Carlo chain in a rather smaller time scale to spot spurious contributions from foreground or other 
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secondary sources. Generalizing the temperature-only estimators we develop estimators for C\ ' for the skew-spectrum (3-point) for specific but 
arbitrary choice of U, V, W, X as well as q( U ' vwx ) an( j q(uv,WX) w j r ^ c j 1 are power spectra of fields related to the specific choice of mixed trispectrum, 
or kurt-spectrum involving temperature and i5-type polarization fields. 

For our next step, we generalized the estimators employed by Yadav, Komatsu,& Wandeltl ilOOH) ; lYadav & Wandeitl J2008h ; lYadav et alj d2008h 
to study the kurt-spectrum. These methods are computationally less expensive and can be implemented using a Monte-Carlo pipeline which involve 
the generation of 3D maps from the cut-sky harmonics using radial integrations of a target theoretical model bispectrum along the line of sight. The 
Monte-Carlo generation of 3D maps is the most computationally expensive part and dominates the calculation. Including polarization field increases 
the computational cost. The technique nevertheless has been used extensively, as it remains highly parallelisable and is near optimal in the presence of 
homogeneous noise and near all-sky coverage. The corrective terms involve linear and quadratic contributions for the lack of spherical symmetry due 
to the presence of inhomogeneous noise and partial sky coverage. These terms can be computed using a Monte-Carlo chain for joint temperature and 
polarization data, but including polarization data requires the ability to handle further complications with an added level of sophistication. We also showed 
that the radial integral involved at the three-point analysis needs to extended to incorporate a double-integral for the mixed trispectrum. For every choice 
of the bi- or trispectrum involving a specific combination of temperature and polarization field, related power spectra can always be defined. These can 
help as a diagnosis for cross-contamination from non-primordial sources in different available frequency channels. 

Though the estimators based on Monte-Carlo analysis based method are very fast, they do not accurately take care of the mode-mode coupling, 
which are present at least at low res olution. We have developed estimators, which are completely optimal even in the presence of inhomogeneous 
no ise and arbitrary sky cover age (e.g. lSmith & Zaldarriagal J2006I) ). These can handle mode-mode c oupling more accurately. Extending previous work 
by iMunshi & Heavens! d2010h which concentrated only on the skew-spectrum, iMunshi et ail d2009h showed how to generalize it to the case of power 
spectrum related to the trispectrum. In this study we have included polarization in a completely general manner both at the level of the bispectrum 
and trispectrum. This involves finding a fast method to construct and invert the joint covariance matrix C; m ;/ m ' in multipole space. In most practical 
circumstances it is possible only to find an approximation to the exact joint covariance matrix, and to cover this we present analysis for an approximate 
matrix which can be used instead of C _1 . This makes the method marginally suboptimal but it remains unbiased. The four-point correlation function also 
takes contributions which are purely Gaussian in nature. The subtraction of these contributions is again simplified by the use of Gaussian Monte-Carlo 
polarized maps with the same power spectrum. A joint Fisher analysis is presented for the construction of the error covariance matrix, allowing joint 
estimation of trispectra contributions from various polarized sources, primaries or secondaries. Such a joint estimation give us fundamental limits on 
how many sources of non-Gaussianity can be jointly estimated from a specific experimental set up which scans the sky for temperature as well as for 
polarization. 

At the level of the bispectrum, primordial non-Gaussianity can for many models be described by a single parameter /jvl- The two degenerate power 
spectra related to the trispectrum we have studied at the 4-point level, require two parameters, typically /jvl and qnl - Use of the two power spectra will 
enable us to put s eparate constraints on /ml and gNL without using information from lower-order analysis of the bispectrum, but they can all be used in 
combination (see Smi dt et all J2010h >. Clearly at even higher-order more parameters will be needed to describe various parameters (/jvl, gNL, Iinl, ■ ■ ■) 
which will all be essential in describing degenerate sets of power spectra associated with multispectra at a specific level. Previous studies concentrated 
only on temperature data, including information from polarisation data can improve the constraints. At present the polarization data is dominated by 
noise, but surveys such as Planck will improve the signal to noise available in polarization data. Future all-sky high sensitivity polarization surveys too 
can further improve the situation and our estimators will be be useful for analysis of such data. Numerical implementation of our estimators will be 
reported elsewhere. We have ignored presence of B-mode polarisation but our formalism can be extended to take into account the magnetic or B-type 
polarisation. 

The analytical results presented here can also be useful in the context of study of shear data from weak lensing surveys. We plan to present our 
results elsewhere (Munshi et al. 2010). 
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